/*
 * fdm1d.cpp
 *
 *  Created on: Aug 29, 2010
 *      Author: nish
 */
#include <iostream>
#include <fstream>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_math.h>

using namespace std;

/* ##### FUNCTION DECLARATIONS #### */
void formHMatrix();
void formNuclPotVector(const double nucLoc);
void formTotalPot();
void printConstants();

/* ################## CONSTANTS AND GLOBAL VARIABLES ############################ */
const double a = 1e-12;
const double L = 1e-9;
const double e=GSL_CONST_MKSA_ELECTRON_CHARGE;
const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
const double m=GSL_CONST_MKSA_MASS_ELECTRON;
const double eps=GSL_CONST_MKSA_VACUUM_PERMITTIVITY;
const double h=GSL_CONST_MKSA_PLANCKS_CONSTANT_H;
const double t0=(hbar*hbar)/(2*m*e); //in eV.
const double estatForce_eps = 1/(4*M_PI*eps);

const double npoints = L/a;

gsl_vector* totalPot;
gsl_vector* nucPot;
gsl_matrix* Hmat;

gsl_vector* eigVal;
gsl_matrix* eigVect;

int main(){


	/*
	 * Calculate potentials
	 */
	//take the nuclues to be present  at the centre -> in the 1D. so, the location will be L/2
	printConstants();
	formNuclPotVector(L/2);

	/*
	 * Form H Matrix
	 */



	/*
	 * Solve for eigenvalues and vectors
	 */

	/*
	 * Create plot files/output files
	 */

	return 0;
}

void printConstants(){
	cout<<"====================================================\n-THE CONSTANTS-\n";
	cout<<"a="<<a<<", L="<<L<<", nucLoc="<<L/2<<" ,Npoints="<<npoints<<endl;
	cout<<"====================================================\n";
}
/*
 *
 */
void formNuclPotVector(const double nucLoc){
	nucPot = gsl_vector_alloc(npoints);
	int i=0;
	double field;
	double dist = 0;
	ofstream outPut;
	char* file = "nucPot.dat";

	gsl_vector_set_zero(nucPot);

	outPut.open(file);

	for(i=0;i<=npoints-1;i++){
		dist += a;
//		if(dist == nucLoc){
//			cout<<"Equals "<<dist<<"\t"<<nucLoc<<endl;
//			field = estatForce_eps*(e*e)/a/1000;
//		}else{
//			field = estatForce_eps*(e*e)/(dist-nucLoc);
//		}
		cout<<i<<", dist="<<fabs(dist-nucLoc)<<" ,field="<<field<<endl;
		field = estatForce_eps*(e)/(fabs(dist-nucLoc));
//		gsl_vector_set(nucPot,i,field);
		outPut<<dist<<"\t"<<field<<endl;
	}
	outPut.close();
}


/**
 *
 */
void formTotalPot(){
	int i=0;
	totalPot = gsl_vector_alloc(npoints);
	double totalPotVal;

	for(i=0;i<npoints;i++){
		totalPotVal = gsl_vector_get(nucPot,i);
		gsl_vector_set(totalPot,i,totalPotVal);
	}
}


/**
 *
 */
void formHMatrix(){
	Hmat = gsl_matrix_alloc(npoints,npoints);
	int i=0;
	gsl_matrix_set_zero(Hmat);

	for(i=0;i<npoints;i++){
		gsl_matrix_set(Hmat,i,i,2*t0 + gsl_vector_get(totalPot,i));
		if(i !=0 && i != npoints-1){
			gsl_matrix_set(Hmat,i,i-1,-t0);
			gsl_matrix_set(Hmat,i,i+1,-t0);
		}else if(i==0){
			gsl_matrix_set(Hmat,i,i+1,-t0);
		}else if(i == npoints-1){
			gsl_matrix_set(Hmat,i,i-1,-t0);
		}
	}
}

